Built with R version 3.6.3
knitr::opts_chunk$set(
echo = TRUE,
message = FALSE,
warning = FALSE,
eval = TRUE,
tidy = TRUE
)
# Environment Set up
rm(list = ls()) #Clean workspace
cat("\014") #Clean Console## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 487947 26.1 1073392 57.4 NA 641482 34.3
## Vcells 926847 7.1 8388608 64.0 32768 1768755 13.5
###################
# Install packages
###################
pkgs <- c("remedy", "Seurat", "dplyr",
"rstudioapi", "cowplot", "ggplot2",
"grid", "gridExtra", "styler", "stringr", "msigdbr")
for(i in 1:length(pkgs)){
if(!require(pkgs[i], character.only = T)){
install.packages(pkgs[i])
require(pkgs[i], character.only = T)
}else{
require(pkgs[i], character.only = T)
}
}
pkgs <- c("gplots", "fgsea", "biomaRt",
"clusterProfiler", "GSEABase", "Bitr",
"org.Mm.eg.db")
for(i in 1:length(pkgs)){
if(!require(pkgs[i], character.only = T)){
BiocManager::install(pkgs[i])
require(pkgs[i], character.only = T)
}else{
require(pkgs[i], character.only = T)
}
}## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'Bitr'
## Warning: package 'Bitr' is not available (for R version 3.6.3)
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'Bitr'
####################
# Github packages
####################
# "EnhancedVolcano" package from Github repo
# devtools::install_github('kevinblighe/EnhancedVolcano')
library(EnhancedVolcano)
# Colour scheme
heatmap.blue.col <- "#2B4FA2"
heatmap.red.col <- "#ED1C24"
white.col <- "#FFFFFF"
DNAM1.neg.col <- "#FF0000" # Red
DNAM1.dim.col <- "#A0A0A4" # Grey
DNAM1.high.col <- "#0000C0" # Blue
col.scheme <- c(DNAM1.neg.col, DNAM1.dim.col, DNAM1.high.col)
# Set working directory to source file location
#setwd(dirname(getActiveDocumentContext()$path))
#working.dir <- getwd()
# create output directories
if(!dir.exists("output")){dir.create("output", recursive = T)}
if(!dir.exists("output/figures")){dir.create("output/figures", recursive = T)}
if(!dir.exists("output/tables")){dir.create("output/tables", recursive = T)}
if(!dir.exists("output/QC")){dir.create("output/QC", recursive = T)}# Load dataset
DNAM1.data <- Read10X(data.dir = "Input_data/10X/")
# Initialize the Seurat object with the raw (non-normalized data). Keep all
# genes expressed in >= 3 cells. Keep all cells with at least 200 detected genes
seurat.object <- CreateSeuratObject(counts = DNAM1.data$`Gene Expression`, min.cells = 3,
min.features = 200)
seurat.object # 13,172 genes and 1,907 cells## An object of class Seurat
## 13172 features across 1907 samples within 1 assay
## Active assay: RNA (13172 features, 0 variable features)
# Get mitochondria genes and calculate percentage of genes
mito.genes <- grep(pattern = "^mt-", x = rownames(x = seurat.object@assays$RNA@data),
value = TRUE)
mito.genes## [1] "mt-Nd1" "mt-Nd2" "mt-Co1" "mt-Co2" "mt-Atp8" "mt-Atp6" "mt-Co3"
## [8] "mt-Nd3" "mt-Nd4l" "mt-Nd4" "mt-Nd5" "mt-Nd6" "mt-Cytb"
percent.mito <- Matrix::colSums(seurat.object@assays$RNA@data[mito.genes, ])/Matrix::colSums(seurat.object@assays$RNA@data)
seurat.object <- AddMetaData(object = seurat.object, metadata = percent.mito, col.name = "percent.mito")
# Visualise mitochondria percentage
VlnPlot(object = seurat.object, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"))## pdf
## 3
## quartz_off_screen
## 2
## 'data.frame': 1907 obs. of 4 variables:
## $ orig.ident : Factor w/ 1 level "SeuratProject": 1 1 1 1 1 1 1 1 1 1 ...
## $ nCount_RNA : num 9729 8060 634 7692 2296 ...
## $ nFeature_RNA: int 2510 2117 411 2230 1045 2140 2589 1564 2767 1849 ...
## $ percent.mito: num 0.0363 0.0362 0.0994 0.0261 0.2134 ...
################### Filter cells
# Remove cells with < 200 genes, and remove cells with > 0.05% mitochondrial
# genes
seurat.object.filt <- subset(seurat.object, subset = nFeature_RNA > 200 & nFeature_RNA <
6000 & percent.mito < 0.05)
seurat.object.filt #13,172 genes and 1,526 cells remain## An object of class Seurat
## 13172 features across 1526 samples within 1 assay
## Active assay: RNA (13172 features, 0 variable features)
# Vis filtering
# Before filtering
VlnPlot(object = seurat.object, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"))# After filtering
VlnPlot(object = seurat.object.filt, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"))## pdf
## 3
## quartz_off_screen
## 2
## An object of class Seurat
## 13172 features across 1907 samples within 1 assay
## Active assay: RNA (13172 features, 0 variable features)
## An object of class Seurat
## 13172 features across 1526 samples within 1 assay
## Active assay: RNA (13172 features, 0 variable features)
# Therefore 381 cells were removed
# Normalise RNA data by logNormalize
seurat.object.filt <- NormalizeData(object = seurat.object.filt, normalization.method = "LogNormalize",
scale.factor = 10000, verbose = TRUE)
# Find variable genes
seurat.object.filt <- FindVariableFeatures(object = seurat.object.filt, verbose = TRUE)
# Scale and regress variables
seurat.object.filt <- ScaleData(object = seurat.object.filt, features = rownames(seurat.object.filt),
vars.to.regress = c("nCount_RNA", "percent.mito"), model.use = "linear", block.size = 2000,
do.scale = TRUE, do.center = TRUE, verbose = TRUE)# Run PCA
seurat.object.filt <- RunPCA(seurat.object.filt, verbose = FALSE)
ElbowPlot(seurat.object.filt, ndims = 50)## pdf
## 3
## quartz_off_screen
## 2
## pdf
## 3
## quartz_off_screen
## 2
seurat.object.filt <- FindNeighbors(seurat.object.filt, dims = 1:20)
seurat.object.filt <- FindClusters(seurat.object.filt, resolution = 0.1)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 1526
## Number of edges: 55479
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9179
## Number of communities: 3
## Elapsed time: 0 seconds
########## tSNE
# tSNE perplex picked based on looped running of RunTSNE at differing values
seurat.object.filt <- RunTSNE(seurat.object.filt, reduction = "pca", dims = 1:20,
perplexity = 30, method = "FIt-SNE", seed.use = 42, dim.embed = 2)
DimPlot(seurat.object.filt, reduction = "tsne") + ggtitle(paste0("Perplexity set at ",
30))## pdf
## 3
## quartz_off_screen
## 2
######### UMAP
# UMAP min.dist and n.neighbors values were picked based on looped running of
# RunUMAP across a range of values
seurat.object.filt <- RunUMAP(object = seurat.object.filt, reduction = "pca", dims = 1:20,
umap.method = "uwot", n.neighbors = 20, min.dist = 0.5, seed.use = 42)
DimPlot(object = seurat.object.filt, reduction = "umap", label = TRUE) + ggtitle(paste0("UMAP Min dist = ",
0.5)) + theme(plot.title = element_text(hjust = 0.5))## pdf
## 3
## quartz_off_screen
## 2
################################## ADD ADT data to seurat object Data needs to be transposed Additionally, cells
################################## filtered during mRNA QC need to be removed
ADT.data <- as.matrix(DNAM1.data$`Antibody Capture`)
dim(ADT.data)## [1] 1 1936
## [1] 1526
# Determine which cells remain in seurat dataset and keep only these from ADT
# dataset
cell.ids.keep <- rownames(seurat.object.filt@meta.data)
keep.logic <- colnames(ADT.data) %in% cell.ids.keep
ADT.data <- ADT.data[, keep.logic]
ADT.data <- t(as.matrix(ADT.data))
rownames(ADT.data) <- "DNAM1"
# Add ADT data to seurat object
seurat.object.filt[["ADT"]] <- CreateAssayObject(counts = ADT.data)
seurat.object.filt## An object of class Seurat
## 13173 features across 1526 samples within 2 assays
## Active assay: RNA (13172 features, 2000 variable features)
## 1 other assay present: ADT
## 3 dimensional reductions calculated: pca, tsne, umap
############################################### Subset cells based on Cd226/DNAM1 ADT expression
# Cells are subsetted by ADT (Cd226/DNAM1) expression into Neg, Dim, High
# populations DNAM1 Neg = ~25% of cells with the lowest expression DNAM1 Dim =
# next ~50% of cells DNAM1 High = top ~25% of cells with highest expression
seurat.object.filt # 1,526 samples/cells## An object of class Seurat
## 13173 features across 1526 samples within 2 assays
## Active assay: RNA (13172 features, 2000 variable features)
## 1 other assay present: ADT
## 3 dimensional reductions calculated: pca, tsne, umap
# Get quantile cutoff values
q.vals <- quantile(t(as.matrix(seurat.object.filt@assays$ADT@data)), probs = c(0.25,
0.75))
print(q.vals)## 25% 75%
## 0.3774276 1.2486700
## [1] 384
sum(seurat.object.filt@assays$ADT@data > q.vals[1] & seurat.object.filt@assays$ADT@data <
q.vals[2]) # 756## [1] 756
## [1] 386
# Get logical variables of which cells fit in each subset
logic.var.neg <- seurat.object.filt@assays$ADT@data <= q.vals[1]
logic.var.dim <- seurat.object.filt@assays$ADT@data > q.vals[1] & seurat.object.filt@assays$ADT@data <
q.vals[2]
logic.var.high <- seurat.object.filt@assays$ADT@data >= q.vals[2]
# Get cell IDs for each subset
DNAM.neg <- colnames(seurat.object.filt@assays$ADT@data)[logic.var.neg]
DNAM.dim <- colnames(seurat.object.filt@assays$ADT@data)[logic.var.dim]
DNAM.high <- colnames(seurat.object.filt@assays$ADT@data)[logic.var.high]
# Create a vector with new subset annotations
cell.ids <- rownames(seurat.object.filt@meta.data)
neg.logic <- cell.ids %in% DNAM.neg
dim.logic <- cell.ids %in% DNAM.dim
high.logic <- cell.ids %in% DNAM.high
DNAM.clust <- cell.ids
DNAM.clust[neg.logic] <- "Neg"
DNAM.clust[dim.logic] <- "Dim"
DNAM.clust[high.logic] <- "High"
DNAM.clust## [1] "Neg" "Dim" "High" "High" "Neg" "Neg" "Dim" "Dim" "Neg" "Neg"
## [11] "Dim" "Dim" "Dim" "High" "Dim" "Dim" "High" "Neg" "Dim" "Dim"
## [21] "Neg" "Dim" "Dim" "Neg" "Neg" "High" "High" "Neg" "Dim" "High"
## [31] "High" "Dim" "Neg" "High" "High" "Dim" "High" "Neg" "Dim" "Neg"
## [41] "Dim" "High" "Neg" "Neg" "Dim" "High" "Neg" "Dim" "Dim" "Dim"
## [51] "Dim" "Dim" "Neg" "Dim" "Dim" "Dim" "Dim" "Neg" "Dim" "Dim"
## [61] "Dim" "Dim" "Dim" "Dim" "Dim" "High" "Dim" "Neg" "Neg" "Neg"
## [71] "Neg" "Dim" "Dim" "High" "Dim" "Dim" "Neg" "High" "High" "High"
## [81] "High" "High" "High" "Dim" "Dim" "Dim" "Neg" "High" "High" "Dim"
## [91] "Dim" "Dim" "Dim" "Dim" "Dim" "Dim" "Neg" "Dim" "High" "High"
## [101] "High" "Dim" "Dim" "High" "Dim" "Dim" "High" "Dim" "Dim" "Dim"
## [111] "Dim" "High" "Dim" "Neg" "Neg" "Dim" "Dim" "High" "Dim" "Dim"
## [121] "High" "Dim" "High" "High" "Dim" "Dim" "Dim" "High" "Neg" "High"
## [131] "Neg" "Dim" "Dim" "Dim" "Dim" "Dim" "Neg" "Dim" "High" "Dim"
## [141] "Dim" "Dim" "High" "Neg" "High" "Dim" "Dim" "Neg" "Neg" "Neg"
## [151] "Dim" "Dim" "High" "Neg" "High" "Neg" "Dim" "Dim" "Neg" "Neg"
## [161] "High" "Dim" "High" "Neg" "Dim" "Neg" "Dim" "Neg" "Dim" "Dim"
## [171] "Dim" "Dim" "Dim" "Dim" "Dim" "Neg" "High" "Dim" "High" "Neg"
## [181] "High" "Neg" "Dim" "Dim" "Neg" "High" "High" "Neg" "High" "Dim"
## [191] "Neg" "Dim" "Dim" "High" "Neg" "High" "Neg" "Neg" "Neg" "Neg"
## [201] "High" "Dim" "Dim" "Neg" "High" "Dim" "Dim" "Dim" "Dim" "High"
## [211] "Neg" "Neg" "Dim" "Dim" "Dim" "Dim" "Dim" "Neg" "Dim" "Dim"
## [221] "Neg" "High" "Dim" "Dim" "Neg" "Dim" "High" "High" "Neg" "High"
## [231] "Dim" "High" "Dim" "High" "Dim" "Neg" "High" "Dim" "High" "Dim"
## [241] "Neg" "Neg" "Dim" "Dim" "Dim" "Dim" "High" "Dim" "Dim" "Dim"
## [251] "High" "Dim" "Dim" "High" "Dim" "Dim" "Neg" "Neg" "Neg" "Dim"
## [261] "High" "Neg" "Dim" "Dim" "Dim" "Dim" "Neg" "High" "Neg" "High"
## [271] "Dim" "Neg" "Neg" "High" "Neg" "Neg" "Dim" "High" "High" "Neg"
## [281] "Neg" "Dim" "High" "Neg" "High" "Dim" "High" "High" "Dim" "Dim"
## [291] "Dim" "High" "High" "Dim" "Dim" "High" "Dim" "High" "High" "Dim"
## [301] "Neg" "Dim" "Dim" "Dim" "High" "Dim" "Neg" "Neg" "Neg" "High"
## [311] "High" "Dim" "Neg" "Neg" "Dim" "Dim" "Dim" "High" "Dim" "Dim"
## [321] "Dim" "Dim" "Neg" "Neg" "Dim" "Dim" "High" "Dim" "Neg" "Dim"
## [331] "High" "Neg" "Dim" "High" "High" "High" "High" "Dim" "Dim" "Neg"
## [341] "Dim" "High" "Dim" "Dim" "High" "High" "Dim" "Neg" "Dim" "Dim"
## [351] "Dim" "Dim" "Dim" "Dim" "Dim" "Dim" "Dim" "Dim" "High" "Neg"
## [361] "High" "Neg" "Dim" "Dim" "Neg" "Neg" "Dim" "Dim" "Dim" "Dim"
## [371] "Dim" "Neg" "Dim" "High" "Dim" "Dim" "Dim" "Dim" "Dim" "Neg"
## [381] "High" "Dim" "Dim" "Dim" "High" "Dim" "High" "Neg" "Neg" "Dim"
## [391] "Neg" "Dim" "Dim" "High" "Dim" "Neg" "Neg" "Dim" "High" "Neg"
## [401] "Neg" "Neg" "Neg" "Dim" "Dim" "Neg" "Neg" "High" "Dim" "Dim"
## [411] "Dim" "Dim" "Neg" "Neg" "Dim" "High" "Dim" "High" "High" "High"
## [421] "Dim" "High" "Dim" "Dim" "Dim" "Neg" "Dim" "High" "Dim" "Neg"
## [431] "High" "Dim" "Dim" "Dim" "Neg" "Dim" "Dim" "Dim" "Dim" "High"
## [441] "Dim" "High" "High" "Neg" "High" "Dim" "High" "Neg" "Dim" "Dim"
## [451] "Neg" "Dim" "Neg" "Dim" "Dim" "Dim" "High" "Neg" "High" "High"
## [461] "Dim" "Neg" "High" "High" "Dim" "Neg" "Dim" "Dim" "Neg" "Dim"
## [471] "Neg" "Dim" "Neg" "High" "Dim" "Neg" "Dim" "Dim" "Dim" "Dim"
## [481] "Dim" "Dim" "Dim" "Dim" "High" "High" "Neg" "High" "Dim" "Dim"
## [491] "Dim" "Dim" "Neg" "Neg" "High" "Dim" "High" "Neg" "Dim" "Dim"
## [501] "Dim" "Neg" "High" "High" "Dim" "Dim" "High" "Dim" "Neg" "Neg"
## [511] "High" "Dim" "Dim" "High" "High" "Dim" "Dim" "Dim" "Dim" "Dim"
## [521] "High" "Neg" "Dim" "Neg" "Neg" "Dim" "Neg" "Dim" "Dim" "High"
## [531] "Neg" "Neg" "Dim" "Dim" "Dim" "Dim" "High" "Dim" "Neg" "High"
## [541] "Dim" "Dim" "High" "High" "Dim" "Dim" "High" "High" "Dim" "Dim"
## [551] "Dim" "High" "Dim" "Dim" "Dim" "High" "Dim" "High" "Dim" "High"
## [561] "Neg" "Dim" "Neg" "Dim" "Dim" "Dim" "Dim" "High" "High" "Dim"
## [571] "Dim" "Dim" "High" "High" "Neg" "High" "Neg" "Dim" "Neg" "Dim"
## [581] "High" "High" "Dim" "Neg" "Dim" "Dim" "High" "High" "Neg" "Neg"
## [591] "Dim" "Dim" "Dim" "Neg" "Dim" "Dim" "Dim" "Dim" "Neg" "Neg"
## [601] "Neg" "Neg" "Dim" "Dim" "Dim" "High" "Dim" "Neg" "Dim" "Dim"
## [611] "High" "Neg" "Dim" "High" "Neg" "Neg" "Neg" "Neg" "Neg" "High"
## [621] "High" "High" "Dim" "Neg" "Neg" "Dim" "Dim" "Dim" "Dim" "Dim"
## [631] "Dim" "Dim" "Neg" "Dim" "Dim" "Dim" "Dim" "Dim" "Dim" "Neg"
## [641] "Neg" "Neg" "Dim" "Dim" "Neg" "High" "High" "Dim" "Neg" "Neg"
## [651] "Dim" "Neg" "High" "Neg" "Dim" "Dim" "Dim" "High" "Dim" "Neg"
## [661] "Dim" "Dim" "Dim" "Dim" "Dim" "Dim" "High" "Dim" "High" "Dim"
## [671] "Neg" "High" "Neg" "Dim" "Dim" "Dim" "Neg" "High" "Neg" "High"
## [681] "Neg" "Neg" "Dim" "Neg" "Dim" "Dim" "Dim" "Dim" "Dim" "High"
## [691] "High" "High" "Dim" "High" "Dim" "High" "Dim" "Neg" "High" "Dim"
## [701] "Dim" "Dim" "Neg" "Dim" "Dim" "High" "Neg" "Dim" "Dim" "High"
## [711] "Neg" "High" "Neg" "Dim" "Dim" "Dim" "Dim" "Dim" "Dim" "Dim"
## [721] "High" "Dim" "Dim" "Dim" "High" "Dim" "High" "Neg" "Dim" "High"
## [731] "Dim" "Neg" "Neg" "Dim" "Dim" "Neg" "Dim" "High" "Dim" "Neg"
## [741] "Dim" "Dim" "Neg" "Dim" "High" "Dim" "High" "Neg" "Dim" "High"
## [751] "Dim" "Dim" "Dim" "Neg" "Dim" "Dim" "Neg" "High" "Dim" "Dim"
## [761] "Dim" "Neg" "Dim" "High" "Dim" "Dim" "Neg" "Neg" "Dim" "Neg"
## [771] "High" "Dim" "Dim" "Dim" "High" "High" "Dim" "Dim" "Dim" "Neg"
## [781] "High" "Dim" "Neg" "Neg" "Dim" "Dim" "Dim" "Neg" "Dim" "High"
## [791] "Dim" "High" "Dim" "High" "High" "Dim" "High" "Dim" "Neg" "Neg"
## [801] "Dim" "Dim" "Dim" "Neg" "Neg" "Neg" "Neg" "High" "Dim" "Neg"
## [811] "Dim" "High" "Dim" "Neg" "Neg" "Dim" "Neg" "Dim" "Dim" "Dim"
## [821] "High" "Neg" "High" "Dim" "Dim" "Neg" "High" "Neg" "Neg" "Dim"
## [831] "High" "Neg" "High" "Dim" "Neg" "Dim" "High" "High" "Neg" "Dim"
## [841] "Dim" "Dim" "Dim" "Dim" "Dim" "Dim" "Dim" "Neg" "Dim" "Neg"
## [851] "Dim" "High" "Neg" "Neg" "Dim" "Dim" "High" "Neg" "High" "High"
## [861] "Dim" "Dim" "Dim" "Neg" "Neg" "Neg" "Dim" "High" "High" "Dim"
## [871] "Dim" "Dim" "Dim" "Dim" "Dim" "Dim" "High" "High" "Dim" "High"
## [881] "Neg" "Neg" "High" "Neg" "Dim" "Neg" "Dim" "High" "High" "High"
## [891] "Neg" "High" "Dim" "Dim" "High" "Dim" "Dim" "Dim" "Neg" "High"
## [901] "High" "Dim" "Dim" "Dim" "Dim" "Neg" "High" "Neg" "High" "Dim"
## [911] "Dim" "High" "Dim" "High" "Dim" "Neg" "Dim" "Dim" "High" "High"
## [921] "High" "High" "High" "High" "High" "High" "High" "High" "Dim" "Dim"
## [931] "Dim" "Dim" "Dim" "Neg" "Dim" "Neg" "Dim" "High" "Dim" "Dim"
## [941] "Dim" "Neg" "Neg" "Dim" "Dim" "Dim" "High" "Dim" "High" "Neg"
## [951] "High" "Dim" "High" "Neg" "High" "Dim" "Neg" "Dim" "Dim" "Dim"
## [961] "Dim" "High" "Dim" "Neg" "Neg" "Dim" "Dim" "Dim" "High" "Dim"
## [971] "Dim" "High" "High" "Dim" "Neg" "Neg" "Dim" "High" "Dim" "Dim"
## [981] "High" "Neg" "High" "High" "Dim" "Neg" "Dim" "High" "Neg" "Neg"
## [991] "Neg" "High" "High" "Dim" "High" "Dim" "High" "Dim" "Neg" "Dim"
## [1001] "High" "Neg" "Neg" "Dim" "Dim" "High" "Neg" "High" "High" "Dim"
## [1011] "Neg" "High" "Dim" "Neg" "Dim" "Dim" "Dim" "Neg" "Dim" "Dim"
## [1021] "Dim" "Dim" "Dim" "High" "Dim" "Dim" "High" "Neg" "High" "Neg"
## [1031] "Neg" "Neg" "Dim" "High" "High" "Neg" "High" "Dim" "Dim" "Neg"
## [1041] "Dim" "Dim" "Dim" "High" "Dim" "High" "Dim" "Neg" "Dim" "Dim"
## [1051] "High" "Dim" "High" "Dim" "Neg" "High" "High" "Dim" "Dim" "Neg"
## [1061] "Neg" "Neg" "Dim" "High" "Dim" "High" "Dim" "Dim" "Dim" "High"
## [1071] "Dim" "Neg" "High" "Neg" "Dim" "Dim" "Dim" "High" "Dim" "Neg"
## [1081] "Neg" "Neg" "High" "Neg" "High" "High" "High" "Dim" "Dim" "Neg"
## [1091] "Dim" "Neg" "High" "Neg" "Dim" "Dim" "Dim" "Neg" "Neg" "Neg"
## [1101] "Neg" "High" "Dim" "Dim" "Dim" "High" "Neg" "Dim" "Neg" "Dim"
## [1111] "Dim" "High" "Neg" "Dim" "High" "Neg" "High" "Neg" "Dim" "Dim"
## [1121] "Dim" "Neg" "High" "Dim" "Dim" "Dim" "High" "Dim" "Dim" "Dim"
## [1131] "High" "Neg" "Dim" "Dim" "Dim" "Neg" "Dim" "Dim" "Neg" "High"
## [1141] "Dim" "High" "Neg" "Dim" "Dim" "Neg" "High" "Dim" "High" "Neg"
## [1151] "Neg" "Neg" "High" "Neg" "Dim" "Dim" "Dim" "Dim" "Dim" "Neg"
## [1161] "Neg" "High" "High" "High" "High" "Dim" "High" "Dim" "Neg" "Dim"
## [1171] "Neg" "Neg" "Dim" "High" "Dim" "Neg" "High" "Dim" "High" "Dim"
## [1181] "High" "High" "Dim" "High" "High" "High" "Neg" "Neg" "High" "High"
## [1191] "Dim" "Dim" "Dim" "Neg" "Dim" "Dim" "Dim" "Dim" "Dim" "Neg"
## [1201] "High" "Neg" "Dim" "Dim" "High" "High" "Dim" "Neg" "Dim" "Dim"
## [1211] "Dim" "Dim" "High" "High" "Neg" "Dim" "Dim" "Dim" "High" "High"
## [1221] "Dim" "Dim" "Dim" "High" "High" "High" "Dim" "Dim" "Dim" "High"
## [1231] "Neg" "Dim" "Neg" "Dim" "Dim" "High" "High" "Dim" "Dim" "Dim"
## [1241] "Neg" "Dim" "Dim" "Dim" "High" "Neg" "Dim" "Dim" "Dim" "Neg"
## [1251] "High" "Dim" "Dim" "High" "High" "Dim" "Dim" "Neg" "Dim" "Neg"
## [1261] "Neg" "High" "Dim" "Neg" "Neg" "Dim" "Dim" "Dim" "High" "High"
## [1271] "Neg" "High" "High" "Dim" "High" "Neg" "Dim" "Dim" "High" "Dim"
## [1281] "High" "Dim" "High" "Dim" "Dim" "Dim" "High" "High" "Dim" "Dim"
## [1291] "Dim" "Dim" "Dim" "High" "Dim" "Neg" "High" "Dim" "Dim" "Neg"
## [1301] "Dim" "Dim" "Dim" "High" "Dim" "High" "Neg" "High" "High" "Dim"
## [1311] "Neg" "Neg" "Dim" "Neg" "Neg" "Neg" "Neg" "Dim" "Dim" "Dim"
## [1321] "Dim" "High" "Dim" "High" "Dim" "High" "Dim" "Dim" "Dim" "Dim"
## [1331] "High" "Dim" "Dim" "Neg" "Dim" "Dim" "High" "Dim" "Dim" "Dim"
## [1341] "High" "Dim" "Dim" "Dim" "Neg" "Neg" "Dim" "High" "Neg" "Dim"
## [1351] "Neg" "High" "High" "High" "Dim" "Neg" "Dim" "Dim" "Neg" "Dim"
## [1361] "Dim" "Dim" "Neg" "Dim" "Neg" "Dim" "Neg" "Dim" "High" "Dim"
## [1371] "Dim" "Dim" "Neg" "Dim" "High" "High" "High" "Neg" "Dim" "Neg"
## [1381] "Dim" "High" "Neg" "Dim" "Dim" "Neg" "Neg" "Neg" "Dim" "Dim"
## [1391] "Dim" "High" "Dim" "Dim" "Neg" "Dim" "Dim" "Neg" "Dim" "Dim"
## [1401] "High" "Dim" "Neg" "Dim" "Dim" "Neg" "Dim" "Dim" "Neg" "Dim"
## [1411] "Dim" "Neg" "Dim" "Dim" "Dim" "Neg" "High" "Dim" "High" "High"
## [1421] "Dim" "Dim" "Dim" "Neg" "Dim" "High" "Neg" "High" "Dim" "High"
## [1431] "High" "Dim" "Dim" "Dim" "Dim" "Neg" "High" "Neg" "High" "Dim"
## [1441] "Dim" "Neg" "Neg" "Dim" "Neg" "Neg" "Neg" "High" "Neg" "Dim"
## [1451] "Dim" "Dim" "High" "Dim" "Dim" "High" "High" "Dim" "Dim" "Dim"
## [1461] "Neg" "High" "Neg" "High" "Neg" "High" "Neg" "Neg" "Dim" "Dim"
## [1471] "High" "Dim" "Dim" "Dim" "Dim" "Neg" "Neg" "Neg" "Dim" "Neg"
## [1481] "Dim" "High" "Dim" "Dim" "Neg" "Neg" "Neg" "Dim" "Neg" "Dim"
## [1491] "Dim" "Dim" "Dim" "Dim" "Neg" "Neg" "Dim" "Dim" "Neg" "Neg"
## [1501] "Dim" "Neg" "High" "Neg" "High" "Dim" "High" "High" "Neg" "Dim"
## [1511] "High" "Dim" "Neg" "Neg" "Neg" "Dim" "Dim" "Neg" "High" "Dim"
## [1521] "Dim" "High" "High" "High" "Neg" "Dim"
# Add ADT clust info to metadata
seurat.object.filt <- AddMetaData(object = seurat.object.filt, metadata = DNAM.clust,
col.name = "DNAM.clust.quantile")
head(seurat.object.filt@meta.data)# Replace ident slot with DNAM clust info
clust.ident <- as.factor(seurat.object.filt@meta.data$DNAM.clust.quantile)
clust.ident## [1] Neg Dim High High Neg Neg Dim Dim Neg Neg Dim Dim Dim High
## [15] Dim Dim High Neg Dim Dim Neg Dim Dim Neg Neg High High Neg
## [29] Dim High High Dim Neg High High Dim High Neg Dim Neg Dim High
## [43] Neg Neg Dim High Neg Dim Dim Dim Dim Dim Neg Dim Dim Dim
## [57] Dim Neg Dim Dim Dim Dim Dim Dim Dim High Dim Neg Neg Neg
## [71] Neg Dim Dim High Dim Dim Neg High High High High High High Dim
## [85] Dim Dim Neg High High Dim Dim Dim Dim Dim Dim Dim Neg Dim
## [99] High High High Dim Dim High Dim Dim High Dim Dim Dim Dim High
## [113] Dim Neg Neg Dim Dim High Dim Dim High Dim High High Dim Dim
## [127] Dim High Neg High Neg Dim Dim Dim Dim Dim Neg Dim High Dim
## [141] Dim Dim High Neg High Dim Dim Neg Neg Neg Dim Dim High Neg
## [155] High Neg Dim Dim Neg Neg High Dim High Neg Dim Neg Dim Neg
## [169] Dim Dim Dim Dim Dim Dim Dim Neg High Dim High Neg High Neg
## [183] Dim Dim Neg High High Neg High Dim Neg Dim Dim High Neg High
## [197] Neg Neg Neg Neg High Dim Dim Neg High Dim Dim Dim Dim High
## [211] Neg Neg Dim Dim Dim Dim Dim Neg Dim Dim Neg High Dim Dim
## [225] Neg Dim High High Neg High Dim High Dim High Dim Neg High Dim
## [239] High Dim Neg Neg Dim Dim Dim Dim High Dim Dim Dim High Dim
## [253] Dim High Dim Dim Neg Neg Neg Dim High Neg Dim Dim Dim Dim
## [267] Neg High Neg High Dim Neg Neg High Neg Neg Dim High High Neg
## [281] Neg Dim High Neg High Dim High High Dim Dim Dim High High Dim
## [295] Dim High Dim High High Dim Neg Dim Dim Dim High Dim Neg Neg
## [309] Neg High High Dim Neg Neg Dim Dim Dim High Dim Dim Dim Dim
## [323] Neg Neg Dim Dim High Dim Neg Dim High Neg Dim High High High
## [337] High Dim Dim Neg Dim High Dim Dim High High Dim Neg Dim Dim
## [351] Dim Dim Dim Dim Dim Dim Dim Dim High Neg High Neg Dim Dim
## [365] Neg Neg Dim Dim Dim Dim Dim Neg Dim High Dim Dim Dim Dim
## [379] Dim Neg High Dim Dim Dim High Dim High Neg Neg Dim Neg Dim
## [393] Dim High Dim Neg Neg Dim High Neg Neg Neg Neg Dim Dim Neg
## [407] Neg High Dim Dim Dim Dim Neg Neg Dim High Dim High High High
## [421] Dim High Dim Dim Dim Neg Dim High Dim Neg High Dim Dim Dim
## [435] Neg Dim Dim Dim Dim High Dim High High Neg High Dim High Neg
## [449] Dim Dim Neg Dim Neg Dim Dim Dim High Neg High High Dim Neg
## [463] High High Dim Neg Dim Dim Neg Dim Neg Dim Neg High Dim Neg
## [477] Dim Dim Dim Dim Dim Dim Dim Dim High High Neg High Dim Dim
## [491] Dim Dim Neg Neg High Dim High Neg Dim Dim Dim Neg High High
## [505] Dim Dim High Dim Neg Neg High Dim Dim High High Dim Dim Dim
## [519] Dim Dim High Neg Dim Neg Neg Dim Neg Dim Dim High Neg Neg
## [533] Dim Dim Dim Dim High Dim Neg High Dim Dim High High Dim Dim
## [547] High High Dim Dim Dim High Dim Dim Dim High Dim High Dim High
## [561] Neg Dim Neg Dim Dim Dim Dim High High Dim Dim Dim High High
## [575] Neg High Neg Dim Neg Dim High High Dim Neg Dim Dim High High
## [589] Neg Neg Dim Dim Dim Neg Dim Dim Dim Dim Neg Neg Neg Neg
## [603] Dim Dim Dim High Dim Neg Dim Dim High Neg Dim High Neg Neg
## [617] Neg Neg Neg High High High Dim Neg Neg Dim Dim Dim Dim Dim
## [631] Dim Dim Neg Dim Dim Dim Dim Dim Dim Neg Neg Neg Dim Dim
## [645] Neg High High Dim Neg Neg Dim Neg High Neg Dim Dim Dim High
## [659] Dim Neg Dim Dim Dim Dim Dim Dim High Dim High Dim Neg High
## [673] Neg Dim Dim Dim Neg High Neg High Neg Neg Dim Neg Dim Dim
## [687] Dim Dim Dim High High High Dim High Dim High Dim Neg High Dim
## [701] Dim Dim Neg Dim Dim High Neg Dim Dim High Neg High Neg Dim
## [715] Dim Dim Dim Dim Dim Dim High Dim Dim Dim High Dim High Neg
## [729] Dim High Dim Neg Neg Dim Dim Neg Dim High Dim Neg Dim Dim
## [743] Neg Dim High Dim High Neg Dim High Dim Dim Dim Neg Dim Dim
## [757] Neg High Dim Dim Dim Neg Dim High Dim Dim Neg Neg Dim Neg
## [771] High Dim Dim Dim High High Dim Dim Dim Neg High Dim Neg Neg
## [785] Dim Dim Dim Neg Dim High Dim High Dim High High Dim High Dim
## [799] Neg Neg Dim Dim Dim Neg Neg Neg Neg High Dim Neg Dim High
## [813] Dim Neg Neg Dim Neg Dim Dim Dim High Neg High Dim Dim Neg
## [827] High Neg Neg Dim High Neg High Dim Neg Dim High High Neg Dim
## [841] Dim Dim Dim Dim Dim Dim Dim Neg Dim Neg Dim High Neg Neg
## [855] Dim Dim High Neg High High Dim Dim Dim Neg Neg Neg Dim High
## [869] High Dim Dim Dim Dim Dim Dim Dim High High Dim High Neg Neg
## [883] High Neg Dim Neg Dim High High High Neg High Dim Dim High Dim
## [897] Dim Dim Neg High High Dim Dim Dim Dim Neg High Neg High Dim
## [911] Dim High Dim High Dim Neg Dim Dim High High High High High High
## [925] High High High High Dim Dim Dim Dim Dim Neg Dim Neg Dim High
## [939] Dim Dim Dim Neg Neg Dim Dim Dim High Dim High Neg High Dim
## [953] High Neg High Dim Neg Dim Dim Dim Dim High Dim Neg Neg Dim
## [967] Dim Dim High Dim Dim High High Dim Neg Neg Dim High Dim Dim
## [981] High Neg High High Dim Neg Dim High Neg Neg Neg High High Dim
## [995] High Dim High Dim Neg Dim High Neg Neg Dim Dim High Neg High
## [1009] High Dim Neg High Dim Neg Dim Dim Dim Neg Dim Dim Dim Dim
## [1023] Dim High Dim Dim High Neg High Neg Neg Neg Dim High High Neg
## [1037] High Dim Dim Neg Dim Dim Dim High Dim High Dim Neg Dim Dim
## [1051] High Dim High Dim Neg High High Dim Dim Neg Neg Neg Dim High
## [1065] Dim High Dim Dim Dim High Dim Neg High Neg Dim Dim Dim High
## [1079] Dim Neg Neg Neg High Neg High High High Dim Dim Neg Dim Neg
## [1093] High Neg Dim Dim Dim Neg Neg Neg Neg High Dim Dim Dim High
## [1107] Neg Dim Neg Dim Dim High Neg Dim High Neg High Neg Dim Dim
## [1121] Dim Neg High Dim Dim Dim High Dim Dim Dim High Neg Dim Dim
## [1135] Dim Neg Dim Dim Neg High Dim High Neg Dim Dim Neg High Dim
## [1149] High Neg Neg Neg High Neg Dim Dim Dim Dim Dim Neg Neg High
## [1163] High High High Dim High Dim Neg Dim Neg Neg Dim High Dim Neg
## [1177] High Dim High Dim High High Dim High High High Neg Neg High High
## [1191] Dim Dim Dim Neg Dim Dim Dim Dim Dim Neg High Neg Dim Dim
## [1205] High High Dim Neg Dim Dim Dim Dim High High Neg Dim Dim Dim
## [1219] High High Dim Dim Dim High High High Dim Dim Dim High Neg Dim
## [1233] Neg Dim Dim High High Dim Dim Dim Neg Dim Dim Dim High Neg
## [1247] Dim Dim Dim Neg High Dim Dim High High Dim Dim Neg Dim Neg
## [1261] Neg High Dim Neg Neg Dim Dim Dim High High Neg High High Dim
## [1275] High Neg Dim Dim High Dim High Dim High Dim Dim Dim High High
## [1289] Dim Dim Dim Dim Dim High Dim Neg High Dim Dim Neg Dim Dim
## [1303] Dim High Dim High Neg High High Dim Neg Neg Dim Neg Neg Neg
## [1317] Neg Dim Dim Dim Dim High Dim High Dim High Dim Dim Dim Dim
## [1331] High Dim Dim Neg Dim Dim High Dim Dim Dim High Dim Dim Dim
## [1345] Neg Neg Dim High Neg Dim Neg High High High Dim Neg Dim Dim
## [1359] Neg Dim Dim Dim Neg Dim Neg Dim Neg Dim High Dim Dim Dim
## [1373] Neg Dim High High High Neg Dim Neg Dim High Neg Dim Dim Neg
## [1387] Neg Neg Dim Dim Dim High Dim Dim Neg Dim Dim Neg Dim Dim
## [1401] High Dim Neg Dim Dim Neg Dim Dim Neg Dim Dim Neg Dim Dim
## [1415] Dim Neg High Dim High High Dim Dim Dim Neg Dim High Neg High
## [1429] Dim High High Dim Dim Dim Dim Neg High Neg High Dim Dim Neg
## [1443] Neg Dim Neg Neg Neg High Neg Dim Dim Dim High Dim Dim High
## [1457] High Dim Dim Dim Neg High Neg High Neg High Neg Neg Dim Dim
## [1471] High Dim Dim Dim Dim Neg Neg Neg Dim Neg Dim High Dim Dim
## [1485] Neg Neg Neg Dim Neg Dim Dim Dim Dim Dim Neg Neg Dim Dim
## [1499] Neg Neg Dim Neg High Neg High Dim High High Neg Dim High Dim
## [1513] Neg Neg Neg Dim Dim Neg High Dim Dim High High High Neg Dim
## Levels: Dim High Neg
seurat.object.filt <- SetIdent(seurat.object.filt, value = clust.ident)
head(seurat.object.filt@active.ident)## AAACCCAGTAATCAAG AAACCCAGTCTGTCCT AAACGCTAGTCACGCC AAAGAACGTTAACAGA
## Neg Dim High High
## AAAGGATAGTCCTACA AAAGGATGTATCGTAC
## Neg Neg
## Levels: Dim High Neg
# Give Active ident and metadata columns factorisation with correct ordering
levels(seurat.object.filt@active.ident)## [1] "Dim" "High" "Neg"
seurat.object.filt@active.ident <- factor(seurat.object.filt@active.ident, levels = c("Neg",
"Dim", "High"))
levels(seurat.object.filt@active.ident)## [1] "Neg" "Dim" "High"
## NULL
seurat.object.filt@meta.data$DNAM.clust.quantile <- factor(seurat.object.filt@meta.data$DNAM.clust.quantile,
levels = c("Neg", "Dim", "High"))
levels(seurat.object.filt@meta.data$DNAM.clust.quantile)## [1] "Neg" "Dim" "High"
##
## Neg Dim High
## 384 756 386
##
## Neg Dim High
## 25.16383 49.54128 25.29489
# Create output directory
if (!dir.exists("output/figures/ADT_subsetting")) {
dir.create("output/figures/ADT_subsetting", recursive = T)
}
# Quick QC
VlnPlot(object = seurat.object.filt, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"))## pdf
## 3
## quartz_off_screen
## 2
# Overall ADT detection Using custom x axis
RidgePlot(seurat.object.filt, "adt_DNAM1", group.by = "orig.ident", log = TRUE, sort = FALSE,
cols = col.scheme) + ggplot2::coord_cartesian(xlim = c(0.9, 4))## pdf
## 3
## quartz_off_screen
## 2
# Vis ADT clustering on VlnPlot
VlnPlot(seurat.object.filt, "adt_DNAM1", group.by = "DNAM.clust.quantile", sort = FALSE,
cols = col.scheme)## pdf
## 3
## quartz_off_screen
## 2
# Vis ADT clustering on RidgePlot
RidgePlot(seurat.object.filt, "adt_DNAM1", group.by = "DNAM.clust.quantile", sort = FALSE,
cols = col.scheme)## pdf
## 3
## quartz_off_screen
## 2
# Visualise Cd226 gene expression
VlnPlot(seurat.object.filt, "Cd226", group.by = "DNAM.clust.quantile", sort = FALSE,
cols = col.scheme)## pdf
## 3
## quartz_off_screen
## 2
RidgePlot(seurat.object.filt, "Cd226", group.by = "DNAM.clust.quantile", sort = FALSE,
cols = col.scheme)## pdf
## 3
## quartz_off_screen
## 2
# DEG of High vs. Neg ADT subsets
DNAMhigh.vs.neg.markers.wilcox <- FindMarkers(seurat.object.filt, ident.1 = "High",
ident.2 = "Neg", logfc.threshold = 0.25, min.pct = 0.1, test.use = "wilcox",
only.pos = FALSE)
write.csv(DNAMhigh.vs.neg.markers.wilcox, "output/tables/DNAMhigh_vs_neg_markers_wilcox.csv")
print(paste0("Number of sig differentially expressed genes = ", sum(DNAMhigh.vs.neg.markers.wilcox$p_val_adj <
0.05)))## [1] "Number of sig differentially expressed genes = 214"
# Do diff gene expression with no threshold or filtering
DNAMhigh.vs.Neg.markers.wilcox.NO.FILT <- FindMarkers(seurat.object.filt, ident.1 = "High",
ident.2 = "Neg", logfc.threshold = 0, min.pct = 0, test.use = "wilcox", only.pos = FALSE)
write.csv(DNAMhigh.vs.Neg.markers.wilcox.NO.FILT, "output/tables/DNAMhigh_vs_Neg_markers_wilcox_ALL_genes_NO_filt.csv")
print(paste0("Number of sig differentially expressed genes = ", sum(DNAMhigh.vs.Neg.markers.wilcox.NO.FILT$p_val_adj <
0.05)))## [1] "Number of sig differentially expressed genes = 900"
save.data.frame.function <- function(df, title) {
x <- as.data.frame(as.matrix(df))
write.table(x, paste0("output/tables/", title, ".txt"), sep = "\t", quote = FALSE)
}
# RNA dataset Write raw data
save.data.frame.function(seurat.object.filt@assays$RNA@counts, "Raw_dataframe_RNA")
# write filtered data
save.data.frame.function(seurat.object.filt@assays$RNA@data, "filtered_dataframe_RNA")
# Write metadata
save.data.frame.function(seurat.object.filt@meta.data, "Meta_data_dataframe")
# ADT dataset Write raw data
save.data.frame.function(seurat.object.filt@assays$ADT@counts, "Raw_dataframe_ADT")
# write filtered data
save.data.frame.function(seurat.object.filt@assays$ADT@data, "filtered_dataframe_ADT")if(!dir.exists("output/figures/VolcanoPlots")){
dir.create("output/figures/VolcanoPlots",
recursive = TRUE)}
#############
# Function
#############
clean.data <- function(x){
# data = data.frame with cols (genes, PValue/FDR, logFC)
colnames(x) <- c("PValue", "logFC", "pct.1", "pct.2", "FDR")
colnames(x)
x$genes <- rownames(x)
head(x)
return(x)
}
# Format data
volcano.data <- clean.data(DNAMhigh.vs.Neg.markers.wilcox.NO.FILT)
# Create colour key for logFC and FDR sig value of interest
FDR.cutoff <- 0.05
logFC.cutoff <- 0.25
# Establish custom colour key for increased vs. decreased genes
# set the base colour as 'black'
keyvals <- rep('grey50', nrow(volcano.data))
# set the base name/label as 'NS'
names(keyvals) <- rep('NS', nrow(volcano.data))
# modify keyvals for variables with fold change > 0.25 & FDR < 0.05
keyvals[which(volcano.data$logFC > logFC.cutoff & volcano.data$FDR < FDR.cutoff)] <- DNAM1.high.col
names(keyvals)[which(volcano.data$logFC > logFC.cutoff & volcano.data$FDR < FDR.cutoff)] <- 'High'
# modify keyvals for variables with fold change < -0.25 & FDR < 0.05
keyvals[which(volcano.data$logFC < -logFC.cutoff & volcano.data$FDR < FDR.cutoff)] <- DNAM1.neg.col
names(keyvals)[which(volcano.data$logFC < -logFC.cutoff & volcano.data$FDR < FDR.cutoff)] <- 'Neg'
unique(names(keyvals))## [1] "High" "Neg" "NS"
# Plot Volcano with top 10 DEGs and custom GOI
# top 10 up and down DEGs
sig.genes <- volcano.data[volcano.data$FDR < 0.05, ]
goi.top <- top_n(sig.genes, 10, logFC)
goi.bottom <- top_n(sig.genes, -10, logFC)
# Hand picked GOIs
curated.goi <- c("Cd226", "Eomes", "Rora",
"Ifng", "Cd69", "Gzma",
"Gzmb", "Cd48", "Gzmk",
"Icam1", "Crtam", "Lag3",
"Tnfrsf9", "Cd247")
curated.goi## [1] "Cd226" "Eomes" "Rora" "Ifng" "Cd69" "Gzma" "Gzmb"
## [8] "Cd48" "Gzmk" "Icam1" "Crtam" "Lag3" "Tnfrsf9" "Cd247"
goi <- unique(c(goi.top$genes, goi.bottom$genes, curated.goi))
# Plot highly customised Volcano
EnhancedVolcano(volcano.data,
lab = rownames(volcano.data),
x = 'logFC',
y = 'FDR',
selectLab = goi,
xlim = c(-1.5, 2.6),
ylim = c(0, 80),
title = "DNAM1 High vs. Neg",
subtitle = "",
ylab = bquote(~-Log[10]~adjusted~italic(P)),
pCutoff = 0.05,
FCcutoff = 0.25,
pointSize = 5,
#col = c("black", "black", "black", "red2"),
colCustom = keyvals,
colAlpha = 0.2,
drawConnectors = TRUE,
widthConnectors = 0.5,
colConnectors = 'grey30',
typeConnectors = "open",
labSize = 4,
gridlines.major = TRUE,
gridlines.minor = FALSE,
border = 'partial',
borderWidth = 0.5,
borderColour = 'black') ## pdf
## 3
## quartz_off_screen
## 2
# Create output directory
if (!dir.exists("output/figures/heatmaps")) {
dir.create("output/figures/heatmaps", recursive = T)
}
# Subset dataset for easier plotting
seurat.object.filt.small <- subset(seurat.object.filt, downsample = 300)
# Top 20 Up- and down-regulated genes in High Vs. Neg comparison
sig.genes <- DNAMhigh.vs.neg.markers.wilcox[DNAMhigh.vs.neg.markers.wilcox$p_val_adj <
0.05, ]
sig.genes$GeneID <- rownames(sig.genes)
up.genes <- top_n(sig.genes, 20, avg_logFC)$GeneID
down.genes <- top_n(sig.genes, -20, avg_logFC)$GeneID
plot.genes <- c(up.genes, down.genes)
# Plot heatmap
DoHeatmap(seurat.object.filt.small, features = unique(plot.genes), assay = "RNA",
group.colors = col.scheme, angle = 90, raster = FALSE)## pdf
## 3
## quartz_off_screen
## 2
## [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
## [6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
## [11] "GO" "GOALL" "IPI" "MGI" "ONTOLOGY"
## [16] "ONTOLOGYALL" "PATH" "PFAM" "PMID" "PROSITE"
## [21] "REFSEQ" "SYMBOL" "UNIGENE" "UNIPROT"
Entrez.ids <- bitr(rownames(DNAMhigh.vs.Neg.markers.wilcox.NO.FILT), fromType = "SYMBOL",
toType = "ENTREZID", OrgDb = "org.Mm.eg.db")
# 6.54% failed to map
# Merge dataset with pval and fold change info
GSEA.dataframe <- merge(Entrez.ids, DNAMhigh.vs.Neg.markers.wilcox.NO.FILT, by.x = "SYMBOL",
by.y = "row.names")
dim(GSEA.dataframe) # 12,127 genes ## [1] 12127 7
# Analysis using msig database
# Msig Database gene sets: H: hallmark gene sets C1: positional gene sets C2:
# curated gene sets C3: motif gene sets C4: computational gene sets C5: GO gene
# sets C6: oncogenic signatures C7: immunologic signatures
# Isolate just the differentially expressed genes
sig.genes <- GSEA.dataframe[GSEA.dataframe$p_val_adj < 0.05, ]
dim(sig.genes) # 869 genes## [1] 869 7
# Subset by upregulated or downregulated genes
up.sig.genes <- sig.genes[sig.genes$avg_logFC > 0, ]
dim(up.sig.genes) # 753 genes## [1] 753 7
# Set variables
cat.val <- "C5"
pval.cut.val <- 0.05
qval.cut.val <- 0.2
# Formatting database
m_t2g <- msigdbr(species = "Mus musculus", category = paste0(cat.val)) %>% dplyr::select(gs_name,
entrez_gene)
# upregulated
e.up <- enricher(up.sig.genes$ENTREZID, TERM2GENE = m_t2g, pvalueCutoff = pval.cut.val,
qvalueCutoff = qval.cut.val)
e.up <- setReadable(e.up, org.Mm.eg.db, keyType = "ENTREZID")
print(head(e.up))## ID
## GO_T_CELL_ACTIVATION GO_T_CELL_ACTIVATION
## GO_LEUKOCYTE_CELL_CELL_ADHESION GO_LEUKOCYTE_CELL_CELL_ADHESION
## GO_LYMPHOCYTE_DIFFERENTIATION GO_LYMPHOCYTE_DIFFERENTIATION
## GO_T_CELL_DIFFERENTIATION GO_T_CELL_DIFFERENTIATION
## GO_REGULATION_OF_T_CELL_ACTIVATION GO_REGULATION_OF_T_CELL_ACTIVATION
## GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION
## Description
## GO_T_CELL_ACTIVATION GO_T_CELL_ACTIVATION
## GO_LEUKOCYTE_CELL_CELL_ADHESION GO_LEUKOCYTE_CELL_CELL_ADHESION
## GO_LYMPHOCYTE_DIFFERENTIATION GO_LYMPHOCYTE_DIFFERENTIATION
## GO_T_CELL_DIFFERENTIATION GO_T_CELL_DIFFERENTIATION
## GO_REGULATION_OF_T_CELL_ACTIVATION GO_REGULATION_OF_T_CELL_ACTIVATION
## GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION
## GeneRatio BgRatio
## GO_T_CELL_ACTIVATION 82/685 448/16881
## GO_LEUKOCYTE_CELL_CELL_ADHESION 61/685 334/16881
## GO_LYMPHOCYTE_DIFFERENTIATION 60/685 340/16881
## GO_T_CELL_DIFFERENTIATION 48/685 242/16881
## GO_REGULATION_OF_T_CELL_ACTIVATION 52/685 315/16881
## GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION 42/685 211/16881
## pvalue
## GO_T_CELL_ACTIVATION 1.854225e-31
## GO_LEUKOCYTE_CELL_CELL_ADHESION 1.677063e-23
## GO_LYMPHOCYTE_DIFFERENTIATION 2.466146e-22
## GO_T_CELL_DIFFERENTIATION 2.865564e-20
## GO_REGULATION_OF_T_CELL_ACTIVATION 3.564098e-18
## GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION 6.108670e-18
## p.adjust
## GO_T_CELL_ACTIVATION 9.211790e-28
## GO_LEUKOCYTE_CELL_CELL_ADHESION 4.165824e-20
## GO_LYMPHOCYTE_DIFFERENTIATION 4.083938e-19
## GO_T_CELL_DIFFERENTIATION 3.559030e-17
## GO_REGULATION_OF_T_CELL_ACTIVATION 3.541288e-15
## GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION 4.650126e-15
## qvalue
## GO_T_CELL_ACTIVATION 7.539865e-28
## GO_LEUKOCYTE_CELL_CELL_ADHESION 3.409733e-20
## GO_LYMPHOCYTE_DIFFERENTIATION 3.342709e-19
## GO_T_CELL_DIFFERENTIATION 2.913072e-17
## GO_REGULATION_OF_T_CELL_ACTIVATION 2.898550e-15
## GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION 3.806136e-15
## geneID
## GO_T_CELL_ACTIVATION Adam8/Anxa1/Apbb1ip/B2m/Bcl3/Camk4/Card11/Casp8/Ccl5/Ccnd3/Ccr2/Cd274/Cd44/Cd47/Cd80/Cd86/Cd8b1/Cdc42/Chd7/Coro1a/Ctla4/Dock2/Dock8/Egr1/Flot2/Gata3/Gpr18/Gpr183/Grap2/Gsn/H2-M3/Havcr2/Icam1/Icos/Ifng/Il18r1/Il2ra/Il7r/Irf1/Irf4/Itgal/Itk/Klrk1/Lck/Lcp1/Lfng/Lgals1/Lgals3/Lgals9/Lilr4b/Map3k8/March7/Msn/Nfkbid/Nfkbiz/Nlrc3/Pag1/Pik3cd/Ppp3ca/Ppp3cb/Prr7/Psmb10/Ptger4/Pycard/Rac2/Rasal3/Rasgrp1/Rhoa/Ripk2/Rora/Runx2/Runx3/Satb1/Sema4a/Smad7/Socs1/Spn/Thy1/Tnfaip8l2/Vsir/Zc3h12a/Zfp683
## GO_LEUKOCYTE_CELL_CELL_ADHESION Adam8/Anxa1/Card11/Ccl5/Ccr2/Cd274/Cd44/Cd47/Cd80/Cd86/Cdc42/Coro1a/Ctla4/Cx3cr1/Dock8/Ezr/Flot2/Gata3/Grap2/H2-M3/Havcr2/Icam1/Icos/Ifng/Il2ra/Il7r/Irf1/Itgal/Itgb1/Itgb2/Itgb7/Klrk1/Lck/Lgals1/Lgals3/Lgals9/Lilr4b/Map3k8/March7/Msn/Nfkbid/Nfkbiz/Pag1/Pycard/Rac2/Rasal3/Rasgrp1/Rhoa/Ripk2/Runx3/Selplg/Skap1/Smad7/Socs1/Spn/Stk10/Thy1/Tnfaip8l2/Vsir/Wnk1/Zc3h12a
## GO_LYMPHOCYTE_DIFFERENTIATION Adam8/Anxa1/B2m/Bcl3/Camk4/Card11/Ccr2/Cd80/Cd86/Chd7/Ctla4/Dock10/Dock2/Egr1/Flt3l/Gata3/Gpr18/Gpr183/H2-M3/Id2/Ifng/Il18r1/Il2ra/Il7r/Inpp5d/Irf1/Irf2bp2/Irf4/Itgb1/Itk/Klf6/Lck/Lfng/Lgals1/Lgals9/Lilr4b/Nfkbid/Nfkbiz/Notch2/Pglyrp2/Pik3cd/Plcl2/Ppp3cb/Prr7/Ptger4/Ptk2b/Rasgrp1/Rhoa/Ripk2/Rora/Runx2/Runx3/Satb1/Sema4a/Smad7/Socs1/Spn/Vsir/Zc3h12a/Zfp683
## GO_T_CELL_DIFFERENTIATION Adam8/Anxa1/B2m/Bcl3/Camk4/Card11/Ccr2/Cd80/Cd86/Chd7/Ctla4/Dock2/Egr1/Gata3/Gpr18/Gpr183/H2-M3/Ifng/Il18r1/Il2ra/Il7r/Irf1/Irf4/Itk/Lck/Lfng/Lgals9/Lilr4b/Nfkbid/Nfkbiz/Pik3cd/Ppp3cb/Prr7/Ptger4/Rasgrp1/Rhoa/Ripk2/Rora/Runx2/Runx3/Satb1/Sema4a/Smad7/Socs1/Spn/Vsir/Zc3h12a/Zfp683
## GO_REGULATION_OF_T_CELL_ACTIVATION Adam8/Anxa1/Camk4/Card11/Ccl5/Ccr2/Cd274/Cd47/Cd80/Cd86/Cdc42/Coro1a/Ctla4/Dock8/Flot2/Gata3/Grap2/Gsn/H2-M3/Havcr2/Icos/Ifng/Il2ra/Il7r/Irf1/Irf4/Klrk1/Lck/Lgals1/Lgals3/Lgals9/Lilr4b/Map3k8/March7/Nfkbid/Nfkbiz/Pag1/Pycard/Rac2/Rasal3/Rasgrp1/Rhoa/Ripk2/Runx3/Smad7/Socs1/Spn/Thy1/Tnfaip8l2/Vsir/Zc3h12a/Zfp683
## GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION Adam8/Anxa1/Card11/Ccl5/Ccr2/Cd274/Cd44/Cd47/Cd80/Cd86/Cdc42/Coro1a/Ctla4/Dock8/Flot2/Gata3/Grap2/H2-M3/Havcr2/Icam1/Icos/Ifng/Il2ra/Il7r/Klrk1/Lck/Lgals1/Lgals9/Lilr4b/Map3k8/Nfkbid/Nfkbiz/Pycard/Rasal3/Rasgrp1/Rhoa/Ripk2/Runx3/Skap1/Socs1/Thy1/Vsir
## Count
## GO_T_CELL_ACTIVATION 82
## GO_LEUKOCYTE_CELL_CELL_ADHESION 61
## GO_LYMPHOCYTE_DIFFERENTIATION 60
## GO_T_CELL_DIFFERENTIATION 48
## GO_REGULATION_OF_T_CELL_ACTIVATION 52
## GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION 42
# Visualisation of enrichment
print(barplot(e.up, showCategory = 10, font.size = 7) + ggtitle(paste0("Upregulated genes msigDB ",
cat.val[i], " enrichment")))## pdf
## 3
## quartz_off_screen
## 2
if (!dir.exists("output/figures/scatterplot")) {
dir.create("output/figures/scatterplot", recursive = T)
}
# Pearson correlation between the two features is displayed above the plot.
FeatureScatter(imputed.seurat, "Cd226", "adt_DNAM1", cols = col.scheme)## pdf
## 3
## quartz_off_screen
## 2
if(!dir.exists("output/figures/VlnPlots_manuscript")){dir.create("output/figures/VlnPlots_manuscript", recursive = T)}
Genes.to.plot <- c("Cd69", "Cd226", "Crtam", "Cd27", "Cd28",
"Tnfrsf4", "Tnfrsf9", "Furin", "Il2ra", "Il12rb2",
"Il7r", "Ccr7", "Cd44", "Mki67", "Cd2",
"Ifng", "Nkg7", "Klrg1", "Slamf1", "Slamf7",
"Klrc1", "S1pr1", "Fasl", "Fas", "Cx3cr1",
"Ctsw", "Cd96", "Entpd1", "Pdcd1", "Adora2a",
"Lag3", "Tgfb1", "Havcr2", "Ctla4", "Cd160",
"Cd274", "Tigit", "Ikzf1", "Ikzf2", "Ikzf3",
"Tox", "Eomes", "Tbx21", "Id2", "Gata3",
"Irf4", "Irf8", "Rora", "Klf2", "Klf3",
"Bach2", "Zeb2", "Myb", "Tcf7", "Itgb2",
"Cdc42", "Itgal", "Itgax", "Itgb1", "Cd48",
"Tln1", "Rap1a", "Rap1b", "Rassf5", "Rac2",
"Ccl5", "Ccl4", "Xcl1", "Il16", "Csf1",
"Gzma", "Gzmb", "Gzmk")
for(i in 1:length(Genes.to.plot)){
print(VlnPlot(imputed.seurat,
paste0(Genes.to.plot[i]),
cols = col.scheme,
#ncol = 3,
pt.size = 0))
dev.copy(pdf, paste0("output/figures/VlnPlots_manuscript/VlnPlot_", Genes.to.plot[i], ".pdf"))
dev.off()
}